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We present a stochastic two-population model that describes the migration and growth of semi- 
sedentary foragers and sedentary farmers along a river valley during the Neolithic transition. The 
main idea of this paper is that random migration and transition from sedentary to foraging way of life 
and backward is strongly coupled with the local crop production and the associated degradation of 
land. We derive a non-linear integral equation for the population density coupled with the equations 
for the density of soil nutrients and crop production. Our model provides an explanation for the 
formation of human settlements along a river valley. The numerical results show that the individual 
farmers have a tendency for aggregation and clustering. We show that the large-scale pattern is a 
transient phenomenon which eventually disappears due to land degradation. 
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I. INTRODUCTION 

The wave of colonization by migrating farmers and establishment of farming communities in Europe 
between 5000 and 3500 BC is currently a topic of great interest in prehistoric archaeology, linguistics and 
anthropology [U, Q • Ammerman and Cavalli-Sforza developed a model for the expansion of farming as a 
demic diffusion which spread into Europe in the form of wave of advance [3, |3| ■ Using the radiocarbon 
dates, they found that farmers spread at an average rate of about one kilometer a year. Interest in 
simulation and spatial modelling of spread of agriculture has been growing rapidly in the last decade, 
especially in the physics community [a, H 0, Q ■ One of the main reasons for this is that the geographical 
spread of population can be effectively described by the classical Fisher-KPP equation and its various 
generalizations @. These models have attracted considerable interest in physics and biology, because 
of the huge number of potential applications. Fort and Mendez have applied a time-delayed theory for 
the Neolithic transition [5] which involves a hyperbolic correction to the Fisher-KPP equation. The 
transition from hunter-gathering to farming did not happen in a uniform way, and that is why Davison 
et al. have taken into account both advection and spatial variation in diffusivity and carrying capacity, 
in the framework of the Fisher-KPP equation Aoki and Shida [ll| have studied the spread of 

farmers into an area occupied by hunter-gatherers in terms of a system of reaction-diffusion equations. 
The main objective of those works was to reproduce the observed rate at which agricultural expansion 
took place in Europe. Despite the interest in establishment of farming communities in Europe, there 
remains no published material on the spatial structure formed behind the wave of advance. The main 
challenge of our work presented below is to set up a model which provides an explanation to not only 
the propagation of farming but also the formation of human settlements. We develop a model that 
demonstrates the tendency of the distribution of population to form clusters, and furthermore, that 
this large-scale pattern in population evolution is a transient phenomenon, which disappears due to 
degradation of land in the form of an extinction wave. Our specific motivation is the successive migration 
of settlements along parallel river valleys in the Tripolye-Cucuteni system, which has been thoroughly 
investigated and documented, e.g. This system can be considered as being one-dimensional, and so 
is particularly amenable to numerical investigation. 



II. A TWO-POPULATION MODEL 



A. A two-population model for semi-sedentary foragers and sedentary farmers 

In our model the population consists of semi-sedentary foragers and sedentary farmers who share the 
same territories. The semi-sedentary foragers are the population of individuals randomly moving from 
place to place along a river valley and searching for food and other resources. An implicit consequence 
of this behaviour is the foundation of new settlements (large localized values of population density) , and 
an interchange between farming and foraging populations. On the contrary the sedentary farmers are 
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individuals who do not migrate. They Hve in smaU villages scattered near cultivated land in the major 
river valleys. Their main activities are the cultivation of soil and crop production. In this paper we 
are interested in the total population density n{x,t) at location x along the river at time t. We define 
this density as n = ni + n2, where rii is the density of semi-sedentary foragers and n2 is the density of 
sedentary farmers. 

We assume that there is not a strict distinction between foraging and sedentary lifestyles. There are 
always random transitions from a sedentary to a foraging way of life, and vice versa; these transitions 
depend strongly on the local food supply. This is one of the main features of our random walk model. 
Regarding the movement of semi-sedentary foragers, they do not jump from place to place completely 
randomly. Unlike Brownian particles in physics, the migration of people cannot be explained by a 
standard diffusion law, in which the flux is proportional to the gradient of number density of individuals. 
To describe the random migration of semi-sedentary foragers and random transitions from one lifestyle 
to another we adopt a biased random walk whose statistical characteristics depend on the local food 
supply. In our model, the probability of a random migration event making a jump z in the time interval 
t to t + At, is XAt. The probability of transition from foraging lifestyle to the farming is aiAt. The 
probability of the conversion of farmers to semi-sedentary foragers is a2At. Thus we introduce a new 
variable, the local crop production per individual per year q{x,t), so that the frequency of jumps A and 
transition rates ai and ck2 depend on the crop production: 

X = X{q), a, ^ ai{q) i = l,2. 

It is natural to assume that the frequency 012 (g) is a decreasing function of q. That is, when sedentary 
farmers are not able to produce enough food to sustain their population, some of them start to migrate 
from their neighborhood at rate X{q{x,t)). 



B. Balance equations for population densities 

We now set up the balance equation for the density of semi-sedentary foragers, ni{x,t). According to 
our random walk model, the density of foragers ni at location x at time t + At can be written as follows 

ni{x,t + At) = {1 - X{q{x,t))At - ai{q{x,t))At)ni{x,t) + 

X{q{x — z, t))Atni{x — z, t)p{z)dz + a2{q{x, t))Atn2{x, t). (1) 

This equation is a conservation law for foragers. The first term on the right hand side represents those 
foragers who stay at location x and do not move during time At and do not become sedentary farmers. 
The second term gives the number of foragers who arrive at x during time At from different places 
X — z, where the jump distance z is distributed by dispersal kernel p (z) . The last term is the number of 
sedentary farmers who become the semi-sedentary foragers during the time At. 

The sedentary farmers do not migrate, and their density n2{x,t) obeys the balance equation involving 
logistic growth and lifestyle transitions: 

n2{x,t + At) = {l-a2{q{x,t))At)n2{x,t) + 

( nix t)\ 

rn2{x,t) (1 - -^J^j At + ai{q{x,t))Atni{x,t), (2) 

where r is the growth rate of the sedentary population. The carrying capacity K is in general, an 
increasing function of the local crop production q. The last term gives the number of foragers who 
become the farmers during the time At. In the limit Ai — s- 0, from ([T]) and ^ we obtain two differential 
equations for ni(x,t) and n2{x,t) : 

Qti f 

-Q^ = X{q{x - z,t))ni{x - z,t)p{z)dz - X{q)ni ~ ai{q)ni + a2{q)n2 (3) 
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aiiq)ni + a2(q)n2. 



(4) 
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In this paper we are interested in the evolution of population density n{x, t) on characteristic time 
scales around 100 — 500 years. Therefore we can adopt a local equilibrium for the two populations, 
describing the balance between the nomadic and sedentary ways of life in proportions p and 1 — p. We 
write 

ni{x,t) = pn{x,t), n2{x,t) ~ {1 — p)n{x,t), (5) 

where the proportion of foragers p is the function of the crop production q : 

p = p{q). (6) 

In the asymptotic regime when a; >> A, one can write the dependence of p on the transition rates ai{q) 
and a2{q) as 

Pil) = — TTT — TT- 
ai[q) + a2(q) 

The evolution equation for biased migration ([3]) and population growth Q can be rewritten as one 
balance equation for the total population density. By adding equations ^ and ^ and using (O we 
obtain a single equation for n(x, t) 



dn 
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p{q{x — z, t))\{q{x — z, t))n{x — z, t)p{z)dz 



-p{q{x, t))\{q{x, t))nix, t) + {1 - p{q))rn ( 1 - ) . (8) 



Now we need to find an approximate form for the function p{q). It is natural to assume that for a large 
crop production q, the population mostly consists of sedentary farmers, that is, p close to zero, while 
for a low crop production the population mostly is semi-sedentary foragers {p is close to one). A simple 
approximation for p might be 

p{q{x, t)) = H {q{x, t) - qmin) , 

where H{x) is the Heaviside step function: 7i(a;) = if a; > 0; H{x) = 1 if a; < 0. Thus farmers do not 
migrate if the crop supply is greater than a critical value q^in- In our numerical simulations we use a 
piece-wise approximation: a linearly decreasing function in the interval between the minimum value qmin 
and the maximum value gmax 

PmSiX: Q — Qmin: 

p{q{x,t)) = { -aq + b, q^in < q < ^max, . 

Pmini Q ^ Qmax 

Based on [l^], we take Qmin, 9max = 300, 736 kg per person per year respectively. For the results discussed 
below we set Pmax = 0.95, Pmin — 0.05 (see Fig. 1) but the exact values are not crucial to our results. 

The dispersal kernel p{z) is assumed to have two cutoffs: a minimum cutoff (foragers do not travel 
very small distances) and a maximum cutoff (foragers do not migrate over very large distances in a single 
jump). In numerical simulations we use minimum cutoff Az of between 5 and 20 km. The maximum 
cutoff is more difficult to estimate. However, since the typical length of the river is around 500 km, it 
seems reasonable to assume that the maximum cutoff is around 50 km (see Fig. 1). 



C. Equation for crop production 

We assume that the human population derives most of their food from the cultivation of land. Thus 
we suggest the following formula for the local food production q{x,t) 

Q{-,t)^a( )(l-e-^^M, (9) 

\no + n[x,tj J \ I 

where F(x,t) is the density of soil nutrients, a is the production rate coefficient, and f3 is the parameter 
that determines how the yield depends on the nutrients. This equation describes how the rate of food 



0.00 -_ 

100 200 300 400 500 o 10 ao 30 40 50 

(a) =! (b) 

Fig. 1: (a) The function q{p) and (b) p{z) for Az = 10 km. 



supply q increases due to the increase in the population density n, and how the degradation of land (the 
decrease of soil nutrients F) leads to a decrease of food production through the factor 1 — e"^-'^^^'*'. Note 
that the factor 1) describes a tendency toward group solidarity that increases the efficiency of food 

production. We know from existing data that the yield, i.e. the production of food per unit area, can 
be up to 0.0736 kg/m^year [H, [3|-We estimate that the area cultivated by people is approximately 10^ 
m^ per person. Thus, we have an estimate for a of 736 kg per person per year. In general, it is difficult 
to model how the production of food depends on the nutrients, since this relation is strongly dependent 
on the environmental conditions. One of the existing models is given by the Mitscherlich-Baule yield 
resp onse [l5j | which takes into account productivity losses due to soil erosion. Following this model, in 
[Iq the factor 1 — e^'^^^^'*^ has been proposed, and this is what we have introduced into our model. The 
value P ~ 890 m^kg^^ is found for the corn cultivation. However, other studies [13] seem to suggest lower 
values (/? ~ 65 m^kg^^). 



D. The land degradation equation 



We adopt the following model for nutrient depletion and the corresponding land degradation. We 
assume that equation for the density of soil nutrients F{x, t) has the form 

^^^^^Cl-^2Fix,t)-jqix,t)n{x,t) (10) 
ot 

where is the rate at which soil nutrients regenerate naturally, ^2 is the rate of nutrient depletion due 
to environmental reasons (erosion, flooding, etc.), 7 is the rate at which nutrients are depleted due to the 
harvests. Here we assume that natural nutrient depletion is much slower than depletion due to human 
activity. In order to estimate 7, we can use the data from [13] ■ This data corresponds to prehistoric 
agriculture in Hawaii; however, similar values are cited for Sub-Saharan agroecosystems and there is no 
data available (as far as we know) for other regions. So using that study we estimate 4 • 10~* kgp/m^yr, 
where kgp denotes kilograms of phosphorus in the soil. This value, together with the data above of 
0.0736 kg/m^yr lead us to assume that the rate of nutrient depletion 7 = 5.43 • 10~^ kgp/kg However, 
other studies mention that the total concentration of nutrients in crops can be up to 3% [l3|- It would 
suggest that the value 7 = 0.03 kgp/kg can be taken as a maximum threshold. Of course, other nutrients 
are also important, but phosphorous can be regarded as a proxy for them. Soil regeneration is a very 
complex process and depends on many environmental parameters. Nevertheless, in some agroeconomics 
models regeneration is modeled as a constant 18] as we have assumed in our model. 



III. NUMERICAL RESULTS 



Numerical simulations of the non-linear integral equation © together with Q and PH)) reveal a 
very interesting dynamical behaviour: the emergence of large-scale patterns in population density. This 
phenomenon can be interpreted as the formation of human settlements along a river valley. Thus our 
model provides an explanation for formation of settlements as a dynamical phenomenon. The individual 
farmers have a tendency for aggregation and clustering as a result of non-linear random migration. 
Moreover our model describes not just a population clustering but subsequent decay of these clusters 
(settlements) due to land degradation. 
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Our model has a rather large number of parameters, but fortunately most of them can be estimated 
from existing sources, as indicated above. For numerical work wc use as unit of time 1 yr, measure 
densities in units of m~^, and assume A = const. We take the values r = 0.03, a = 736, /3 = 200, 
K = 10~^, ^1 = 10~^, ^2 = 0, qmin = 300, QmaK = 736. Note that for simplicity in this initial study 
we take iiT to be a constant, rather than a function of q. We explore ranges 0.01 < 7 < 0.03 yr""'^, 
0.05 < A < 0.3 yr^^, 5 < Az < 20 km. Our initial conditions were rather arbitrary: n = O.IK in 
< a; < 15 km, n = elsewhere, although from experiments with other initial conditions with n(x, 0) 
non-zero in a: < 100 it appears that results are quite insensitive to details of the initial state. Again 
for simplicity, we set the initial value F{x,0) to be uniform and equal to Fq. Our 'standard' choice was 
Fo = 0.01. If the range of x were to be extended to negative values, the evolution of n{x,t) would be 
approximately symmetric about x = 0. 

Fig. 2 shows the population density n{x,t) for A = 0.1, with Ax = 10, 7 = 0.01, at intervals of 100 
yr. Panel (a) shows the "arrival" of farmers at location < .t < 15 along the river. After 300 years, 
one can see clear evidence of clustering of population (" settlements" ) . Panel (e) shows their subsequent 
decay and reappearance of clusters behind the extinction wave. The entire population decays over about 
800 yr, depending to some extent on the exact choice of parameters. A general feature of our modelling 
is that most population clusters grow and decay in sit,u, without significant movement. We can make an 
order of magnitude estimate of the total population of a cluster by taking the linear extent of a cluster 
as the diameter of a circular settlement. This gives typical figures of 0(10"^) individuals. 

We also explored the dependence of the speed of population advance, and the separation of the clusters, 
as a function of the parameters. The speed of advance is rather ill-defined, but for simplicity we monitored 
the movement of the main peak of the distribution. (Arguably we could have, e.g., studied the position of 
the cluster furthest from the origin.) We found that the speed of advance depends approximately linearly 
on both A and 7, and to be insensitive to Ax. The separation of clusters is independent of both A and 
7 and depends approximately linearly on Ax (= (1.5 — 2)Ax). If A < 0.05 or Ax > 15 the clustering 
phenomenon does not occur. While we have concentrated our eflforts on studying variation among a 
manageable subset of parameters while retaining fixed plausible values for the others, we did also make 
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a small study of the effects of varying Fq and ^i. Increasing the value of Fq, up to 0.1 from our canonical 
value of 0.01, prolongs the timescale of the population evolution, by a factor of up to about 2. Values of 
Fq significantly smaller than 0.01 do not give clustering for typical choices of the other parameters. It 
means that for non-fertile soils transition to farming and the emergence of clusters are unlikely to occur. 
For larger values of ^i, > 10~^, the regeneration is so strong that, although clustering occurs initially 
much as for = lO""*, the final state is a spatially uniform population. For smaller values of ^i, e.g. 
^1 = 10~^, we still see strong clustering, but the overall phenomenon persists for a significantly shorter 
time, as the episodes of cluster regeneration seen when = 10~* do not now occur. 

IV. DISCUSSION AND CONCLUSIONS 

We have explored a continuous model of linear (one-dimensional) population migration and clustering. 
We took our inspiration and guidance from the successive south to north migration of the Tripolye- 
Cucuteni cultures along parallel river valleys ^l5|. We attempted only to illustrate the migration along 
a single valley. The main result is that our model gives an explanation for the formation of settlements 
as a dynamical phenomenon. The individuals have a tendency for aggregation and clustering as a re- 
sult of non-linear random migration. Moreover our model describes subsequent decay of these clusters 
(settlements) due to land degradation in the form of an extinction wave. Clearly our model, being math- 
ematically continuous, cannot be expected to reproduce in detail the essentially discrete phenomenon of 
the establishment, growth and decay of major settlements and their satellites, as revealed by the archae- 
ological record - our goal is far more modest. For reasonable choices of parameters we obtain migration 
over distances of order 500 km in times of 500 — 1000 yr. Without fine tuning, our solutions exhibit 
distinct clustering of population, at intervals of 10 — 30 km. Plausibly these clusters are estimated to 
have O(IO^) individuals. These estimates are all consistent with what is known about the Tripolye- 
Cucuteni cultures [13|]. Overall, we feel that our quite naive modelling captures some of the essence of 
the clustering seen in population migration, and points the way for more sophisticated modelling. We 
can think of a number of significant developments in the future, including taking K — K(q) and allowing 
for a latency effect, in which population in a cluster ('settlement') has a reduced probability of moving, 
representing an attachment to the investment of building a house and clearing land. It is easy to allow 
for additional resources available from e.g. forest or river, which are maybe harder to exhaust than the 
fertility of the land and so can act as a reservoir of resource. It would be interesting to apply our model 
to understand societal collapses to which environmental problems as habitat destruction, soil degradation 
and overpopulation contribute p^ . 
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